library(dplyr)
library(stringr)
library(ggplot2)
library(scales)
library(ggrepel)
library(ggpubr)

if (!file.exists("output")) {
  dir.create("output")
}

# Bubble plot displaying associations between tobacco-related signatures and anatomical subsites

# Signatures and variables of interest
signatures <- c("SBS4", "SBS92", "SBS_I", "DBS2", "DBS6", "ID3", "DBS_D", "ID4")
vars <- c("Larynx yes", "OC yes", "OPC yes", "Hypopharynx yes") %>% rev()

# NOTE: Output obtained from risk_factor_regressions.R, also shown in manuscript Supplementary Table 13
bubbleplot.tob <- read.csv("ExtendedDataFigure_5ab_data.csv") %>%
  mutate(
    p_adj_log = -log10(p_adj),
    sig_type = case_when(
      grepl("SBS", signature) ~ "SBS",
      grepl("DBS", signature) ~ "DBS",
      grepl("ID", signature) ~ "ID"
    )
  ) %>%
  filter(
    Independent_variable %in% vars,
    signature %in% signatures
  ) %>%
  mutate(
    signature = factor(signature, levels = signatures),
    log_OR = log2(OR),
    signif = ifelse(p_adj < 0.05, "p adj < 0.05", "p adj > 0.05")
  )

ggplot(bubbleplot.tob, aes(
  y = Independent_variable, x = signature,
  color = log_OR, size = p_adj_log,
  shape = signif
)) +
  geom_point(stroke = 0.75) +
  scale_x_discrete(limits = signatures) +
  scale_y_discrete(limits = vars) +
  scale_color_gradient2(high = "#ca0020", mid = "white", low = "#0571b0", midpoint = 0, limits = c(-4, 4), breaks = c(-4, -2, 0, 2, 4)) +
  scale_size_continuous(breaks = c(0, 3, 6, 9), limits = c(0, 9)) +
  scale_shape_manual(values = c(19, 21)) +
  labs(
    title = "Associations with subsite",
    subtitle = "Corrected for sex, age of diagnosis, region, tobacco, \nand alcohol exposure",
    x = "", y = "", color = "log2(OR)", size = "-log10(p adj)", shape = ""
  ) +
  theme_bw() +
  theme(
    legend.position = "right",
    legend.text = element_text(size = 12),
    plot.title = element_text(face = "bold", size = 12),
    plot.subtitle = element_text(size = 9),
    axis.text = element_text(size = 12),
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank()
  )

ggsave(last_plot(),
  filename = paste0("output/ExtendedData_5c_subsites_", Sys.Date(), ".pdf"),
  device = "pdf",
  width = 6.5, height = 3.5, dpi = 700
)
